function cfdAssembleDiffusionTerm(theEquationName)
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018 
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   This function assembles Diffusion term for scalar equation
%--------------------------------------------------------------------------

% Assemble over Interior Faces
cfdAssembleDiffusionTermInterior(theEquationName);

% Assemble vver Boundary Patches
theNumberOfBPatches = cfdGetNumberOfBPatches;
for iBPatch=1:theNumberOfBPatches

    % Find the Physical Type
    theBoundary = cfdGetBoundaryPatchRef(iBPatch);
    thePhysicalType = theBoundary.type;
    theBCType = cfdBcForBoundaryPatch(iBPatch);

    if strcmp(thePhysicalType,'wall')
        if strcmp(theBCType,'zeroGradient')
            cfdAssembleDiffusionTermZeroGradient(theEquationName, iBPatch);
        elseif strcmp(theBCType,'fixedValue')
            cfdAssembleDiffusionTermSpecifiedValue(theEquationName, iBPatch);
        else
            error([theBCType ' bc is not implemented']);
        end
    elseif strcmp(thePhysicalType,'inlet')
        if strcmp(theBCType,'fixedValue')
            cfdAssembleDiffusionTermSpecifiedValue(theEquationName, iBPatch);
        elseif strcmp(theBCType,'zeroGradient')
            cfdAssembleDiffusionTermZeroGradient(theEquationName, iBPatch);
        else
            error([theBCType ' bc is not implemented']);
        end
    elseif strcmp(thePhysicalType,'outlet')
        if strcmp(theBCType,'fixedValue')
            cfdAssembleDiffusionTermSpecifiedValue(theEquationName, iBPatch);
        elseif strcmp(theBCType,'zeroGradient')
            cfdAssembleDiffusionTermZeroGradient(theEquationName, iBPatch);
        else
            error([theBCType ' bc is not implemented']);
        end
    elseif strcmp(thePhysicalType,'symmetry') || strcmp(thePhysicalType,'empty')
        continue;
    else
        error([thePhysicalType ' physical condition is not implemented']);
    end    
end

end




%===================================================
% INTERIOR
%===================================================
function cfdAssembleDiffusionTermInterior(theEquationName)

% Get info
theNumberOfInteriorFaces = cfdGetNumberOfInteriorFaces;
Sf = cfdGetFaceSfSubArrayForInterior;% 内部面的surface vector
CF = cfdGetFaceCFSubArrayForInterior;% 内部面的从owner指向neighbor的向量
owners_f = cfdGetOwnersSubArrayForInteriorFaces;% 内部面的owner索引序列
neighbours_f = cfdGetNeighboursSubArrayForInteriorFaces;% 内部面的neighbor索引序列
iFaces = 1:theNumberOfInteriorFaces;% 内部面索引序列

% Calculated info
e = cfdUnit(CF);% 内部面的从owner指向neighbor的单位化向量
cfdMagSf = cfdMag(Sf);
cfdMagCF = cfdMag(CF);

% Get fields 
phi = cfdGetSubArrayForInterior(theEquationName);% 得到内部element处的物理量的场
gradPhi = cfdGetGradientSubArrayForInterior(theEquationName);% 得到内部element处gradient的场
gradPhi_f = cfdInterpolateGradientsFromElementsToInteriorFaces('linear',gradPhi);% 将位于element位置的梯度线性插值到内部面上

% Get gamma
gamma = cfdGetSubArrayForInterior('gamma');% 得到内部element处的扩散系数
gamma_f = cfdInterpolateFromElementsToInteriorFaces('linear',gamma);% 将定义于element位置的扩散系数插值到内部面上

% Minimum correction approach（但实际上看着是Orthogonal Correction的方法？？？）
Ef = [cfdMagSf.*e(:,1),cfdMagSf.*e(:,2),cfdMagSf.*e(:,3)];

% Calculate non-orthogonal complement of Sf
Tf = Sf - Ef;

% Geometric diffusion
geoDiff_f = cfdMag(Ef)./cfdMagCF;

% Linear fluxes
local_FluxCf =   gamma_f.*geoDiff_f;
local_FluxFf = - gamma_f.*geoDiff_f;

% Non-linear fluxes
local_FluxVf = - gamma_f.*dot(gradPhi_f(:,:)',Tf(:,:)')';% 将内部面上的梯度点乘非正交部分Tf，纳入源项计算

% Update global fluxes
theFluxes = cfdGetFluxes;

theFluxes.FluxCf(iFaces,1) = local_FluxCf;
theFluxes.FluxFf(iFaces,1) = local_FluxFf;
theFluxes.FluxVf(iFaces,1) = local_FluxVf;
% FluxTf是整个这个面的通量
theFluxes.FluxTf(iFaces,1) = theFluxes.FluxCf(iFaces).*phi(owners_f) + theFluxes.FluxFf(iFaces).*phi(neighbours_f) + theFluxes.FluxVf(iFaces);

cfdSetFluxes(theFluxes);

end

%===================================================
% Specified Value
%===================================================
function cfdAssembleDiffusionTermSpecifiedValue(theEquationName, iBPatch)

% Get info
Sf_b = cfdGetFaceSfSubArrayForBoundaryPatch(iBPatch);
wallDist_b = cfdGetWallDistSubArrayForBoundaryPatch(iBPatch);
owners_b = cfdGetOwnersSubArrayForBoundaryPatch(iBPatch);

iBFaces = cfdGetBFaceIndicesForBoundaryPatch(iBPatch);

% Get required fields
gamma_b = cfdGetSubArrayForBoundaryPatch('gamma', iBPatch);

phi = cfdGetDataArray(theEquationName);
phi_C = phi(owners_b,:);
phi_b = cfdGetSubArrayForBoundaryPatch(theEquationName, iBPatch);

% Geometric diffusion 直接令Tb=0，Sb=Eb？这跟书上式(8.81)(8.82)那样拆分出Eb和Tb不同
geoDiff_b = cfdMag(Sf_b)./wallDist_b;

% Linear fluxes 正交项
local_FluxCb =  gamma_b.*geoDiff_b;
local_FluxFb = -gamma_b.*geoDiff_b;

% non-linear fluxes 非正交项
local_FluxVb = zeros(size(local_FluxFb));

% Update global fluxes
theFluxes = cfdGetFluxes;

theFluxes.FluxCf(iBFaces) = local_FluxCb;
theFluxes.FluxFf(iBFaces) = local_FluxFb;
theFluxes.FluxVf(iBFaces) = local_FluxVb;
theFluxes.FluxTf(iBFaces) = local_FluxCb.*phi_C + local_FluxFb.*phi_b + local_FluxVb;

cfdSetFluxes(theFluxes);

end


%===================================================
% Zero Gradient
%===================================================
function cfdAssembleDiffusionTermZeroGradient(theEquationName, iBPatch)

end